clear
*Analyse the welfare simulation results and show the most important indicator results: 

frame create three
frame create three_elast
frame create fee_schedule


frame change three
use "$root/Results/estimation_results/welfare_max08_3y_1.dta", clear

forvalues i=2/6 {
append using "$root/Results/estimation_results/welfare_max08_3y_`i'.dta"
}


*Get the initial levels to be able to calculate the differences: 
sort tax_sim subsidy 
gen CO2_d=(co2_ann-co2_ann[1])/1000
gen emb_CO2_d=(emb_co2_ann-emb_co2_ann[1])/1000
gen tot_CO2_d=CO2_d+emb_CO2_d
gen petrol_tax_d=petrol_tax-petrol_tax[1]
gen car_tax_d=tax_total-tax_total[1]
gen tax_pv_d=tax_pv-tax_pv[1]
gen veh_tax_d=veh_tax-veh_tax[1]

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen pub_d=(-paid_subsidy+veh_tax_d+tax_pv_d+(petrol_tax_d)*`cap')/1000

gen wel_d=(CSexp_change*1000-paid_subsidy+veh_tax_d+tax_pv_d+(petrol_tax_d+abs(tot_CO2_d)*175)*`cap')/1000
replace wel_d=0 if missing(wel_d)

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen co2_pv=CO2_d*`cap'

gen abat_costs=1000*(CSexp_change+pub_d)/co2_pv

gen mon_abat_costs=pub_d/co2_pv

frame change fee_schedule
use "$root/Data/Original/tax_params_ext.dta", clear

frame change three
sort tax_sim subsidy
gen CO2_perc=(co2_ann-co2_ann[1])*100/co2_ann[1]


frlink m:1 tax_sim, frame(fee_schedule)

frget _all, from(fee_schedule)

**Generate the outcome which is the absolut change in public budget: 

gen abs_pub_d=abs(pub_d)

**Check the outcome with cost deviation minimization and experienced consumer surplus
quietly su abs_pub_d if CSexp_change>=0 & prob_EV>=0.0231
list subsidy tax_sim prob_EV CS_change be_change CSexp_change CO2_d CO2_perc pub_d car_tax_d paid_subsidy *_K if abs_pub_d==`r(min)'

**Check the outcome with cost deviation minimization and weighted experienced consumer surplus
quietly su abs_pub_d if weight_CSexp_d>=0 & prob_EV>=0.0231
list subsidy tax_sim prob_EV CS_change be_change CO2_d CO2_perc pub_d car_tax_d paid_subsidy *_K if abs_pub_d==float(`r(min)')

**What if we weighted with squared wealth?
quietly su abs_pub_d if weight_CSexp_d_sq>=0 & prob_EV>=0.0231
list subsidy tax_sim prob_EV CS_change be_change CO2_d CO2_perc pub_d car_tax_d paid_subsidy *_K if abs_pub_d==float(`r(min)')


*Do it for the elastic one: 
frame change three_elast
use "$root/Results/estimation_results/welfare_max08_3y_elast_1.dta", clear

forvalues i=2/7 {
append using "$root/Results/estimation_results/welfare_max08_3y_elast_`i'.dta"
}

gen elast=-0.3


*Get the initial levels to be able to calculate the differences: 
sort tax_sim subsidy 
gen CO2_d=(co2_ann-co2_ann[1])/1000
gen emb_CO2_d=(emb_co2_ann-emb_co2_ann[1])/1000
gen tot_CO2_d=CO2_d+emb_CO2_d
gen petrol_tax_d=petrol_tax-petrol_tax[1]
gen car_tax_d=tax_total-tax_total[1]
gen tax_pv_d=tax_pv-tax_pv[1]
gen veh_tax_d=veh_tax-veh_tax[1]

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen pub_d=(-paid_subsidy+veh_tax_d+tax_pv_d+(petrol_tax_d)*`cap')/1000

gen wel_d=(CSexp_change*1000-paid_subsidy+veh_tax_d+tax_pv_d+(petrol_tax_d+abs(tot_CO2_d)*175)*`cap')/1000
replace wel_d=0 if missing(wel_d)

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen co2_pv=CO2_d*`cap'

gen abat_costs=1000*(CSexp_change+pub_d)/co2_pv

gen mon_abat_costs=pub_d/co2_pv


frlink m:1 tax_sim, frame(fee_schedule)

frget _all, from(fee_schedule)

**Generate the outcome which is the absolut change in public budget: 

gen abs_pub_d=abs(pub_d)

**Check the outcome with cost deviation minimization and experienced consumer surplus
quietly su abs_pub_d if CSexp_change>=0 & prob_EV>=0.0231
list subsidy tax_sim prob_EV CSexp_change CS_change be_change CO2_d emb_CO2_d tot_CO2_d pub_d car_tax_d paid_subsidy *_K if abs_pub_d==float(`r(min)')

**Check the outcome with cost deviation minimization and weighted experienced consumer surplus
quietly su abs_pub_d if weight_CSexp_d>=0 & prob_EV>=0.0231
list subsidy tax_sim prob_EV CSexp_change CS_change be_change tot_CO2_d pub_d car_tax_d paid_subsidy *_K if abs_pub_d==float(`r(min)')

**What if we weighted with squared wealth?
quietly su abs_pub_d if weight_CSexp_d_sq>=0 & prob_EV>=0.0231
list subsidy tax_sim prob_EV CSexp_change CS_change be_change tot_CO2_d pub_d car_tax_d paid_subsidy *_K if abs_pub_d==float(`r(min)')


**It mostly stays the same - We now recalculate the "optimal" scenarios to be able to also present the results for each weatlh group:


frame change default


*read in the data: 
use "$root/Data/Produced/EV_welfare04.dta", clear

order price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4, after(choice)

replace var_cost=0 if missing(var_cost)

sort id car_option

*Set up the matrices and vectors needed for the computation:
mata

ndraws=150

b_f= (0.021183645,-2.593121925,0.154611204,-0.321820168,0.037136405,0.167413791,0.378321531,0.614503527,0.858999103,-0.005445116,-0.010654642,0.136554812,-0.083769393,-0.022510516,0.007285076,0.06246904,1.476449254,2.260850432,0.424973893,0.853075696,0.6579357,1.4180835,0.027265514,0.980066634,0.090419384,0.171157248,-0.142510175,-0.661616614,-0.177476937,0.595728981,0.214934965,0.378106038,-0.080343758,0.394597341,0.734410779,0.078025015,1.019999496,0.291640034,-0.620171333,-0.126470184,-0.114005609,-0.169462474,0.335901857,-0.232550483,0.070071269,0.05306741,0.003621919,0.006019563,0.024778481)'
/*
** Define the estimated vectors of random coefficients:
*/
beta_f=J(1,ndraws,b_f)
rseed(1234)
beta_p = rnormal(1,ndraws,-0.080965469,0.013364931)
beta_he = rnormal(1,ndraws,0.263640717,0.000661958)
beta_we = rnormal(1,ndraws,0.000192522,0.000045627)
beta_var = rnormal(1,ndraws,-0.040212482,0.001878464)


beta= (beta_p \ beta_he \ beta_we \ beta_var \ beta_f)

mata drop beta_he beta_we b_f


wealth_group=.
st_view(wealth_group,.,"wealth_quart2 wealth_quart3 wealth_quart4")


/*
Define the price coefficients:
*/
b_p=(0.003621919,0.006019563,0.024778481)'


b_w=wealth_group*b_p

mata drop wealth_group b_p



b_price=J(rows(b_w),1,beta_p) + J(1,ndraws,b_w)

gamma=beta_var:/b_price

mata drop beta_p b_w beta_var

end



**To be able to calculate the changes in CS we need the initial levels: 
frame create cs_init
frame change cs_init
import excel "$root/Results/estimation_results/counterfactuals/CS3.xlsx", firstrow

frame change default


**
*First case: Elasticity=0, tax_sim=12, Subsidy=6200
**
local count=1
local y=0
local x=12
local i=6200


**This is the program to calculate the optimal outcomes by quartile: 

gen elast=`y'

gen car_tax_alt=car_tax
replace car_tax_alt=car_tax/0.4 if drivetype==3
replace car_tax_alt=car_tax/0.6 if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax/0.8 if drivetype!=3 & env_cat==2	
	
	gen tax_sim=`x'
	
	merge m:1 tax_sim using "$root/Data/Original/tax_params_ext.dta", 
	drop if _merge==2
	drop _merge
	
	
local r=0.06

gen pv_tax_alt=pv_tax

*for EV
replace pv_tax_alt=car_tax_alt*(1+EV_K)+(car_tax_alt*(1+EV_K))/((1+`r')^(1))+(car_tax_alt*(1+EV_K))/((1+`r')^(2))+(car_tax_alt*(1+EV_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype==3 

*for category A
replace pv_tax_alt=car_tax_alt*(1+A_K)+(car_tax_alt*(1+A_K))/((1+`r')^(1))+(car_tax_alt*(1+A_K))/((1+`r')^(2))+(car_tax_alt*(1+A_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

*for category B
replace pv_tax_alt=car_tax_alt*(1+B_K)+(car_tax_alt*(1+B_K))/((1+`r')^(1))+(car_tax_alt*(1+B_K))/((1+`r')^(2))+(car_tax_alt*(1+B_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*for category E
replace pv_tax_alt=car_tax_alt*(1+E_K)+(car_tax_alt*(1+E_K))/((1+`r')^(1))+(car_tax_alt*(1+E_K))/((1+`r')^(2))+(car_tax_alt*(1+E_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==5
*for category F
replace pv_tax_alt=car_tax_alt*(1+F_K)+(car_tax_alt*(1+F_K))/((1+`r')^(1))+(car_tax_alt*(1+F_K))/((1+`r')^(2))+(car_tax_alt*(1+F_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax_alt*(1+G_K)+(car_tax_alt*(1+G_K))/((1+`r')^(1))+(car_tax_alt*(1+G_K))/((1+`r')^(2))+(car_tax_alt*(1+G_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==7

replace car_tax_alt=car_tax_alt*(1+G_K) if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax_alt*(1+F_K) if drivetype!=3 & env_cat==6
replace car_tax_alt=car_tax_alt*(1+E_K) if drivetype!=3 & env_cat==5
replace car_tax_alt=car_tax_alt*(1+B_K) if drivetype!=3 & env_cat==2
replace car_tax_alt=car_tax_alt*(1+A_K) if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax_alt*(1+EV_K) if drivetype==3 


gen ann_cost=drive_cost_ann+car_tax
gen ann_cost_alt=drive_cost_ann+car_tax_alt


gen alt_km=average_use*(1+elast*((ann_cost_alt-ann_cost)/ann_cost))

replace alt_km=average_use if missing(alt_km)

gen alt_cost_drive=(alt_km/100)*drive_cost
gen alt_drive_cost=alt_cost_drive+alt_cost_drive*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost4=(alt_drive_cost+pv_tax_alt)/1000

replace alt_var_cost4=0 if missing(alt_var_cost4)


drop alt_drive_cost alt_cost_drive ann_cost ann_cost_alt

	
		gen subsidy=`i'
	
	gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-subsidy if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

drop price_alt_tot

sort id car_option

*Add the mata program here: 

mata {
X=. ;
st_view(X,.,"price_alt height weight alt_var_cost4 typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4") ;

var_cost=.
st_view(var_cost,.,"alt_var_cost4")

l_fit=exp(X*beta) ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
fit=mean(l_fit')
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

/*
mata drop prob_all be var_cost 
*/
prob_all=. ;
be=. ;
var_cost=. ; 

st_addvar("double", "prob_alt")
st_store(.,"prob_alt",prob')

st_addvar("double", "be_alt")
st_store(.,"be_alt",be_m')


CS=mean((((-1):/b_price):*log(bigtot))')

/*
mata drop bigtot
*/
bigtot=. ; 

st_addvar("double", "CS_alt")
st_store(.,"CS_alt",CS')

/*
mata drop prob CS be_m
*/
prob=.;
CS=.;
be_m=.;

}


frame put prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt drivetype subsidy prhnewprice tax_sim elast car_option wealth_group, into(probs)

frame change probs

gcollapse (mean) prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt prhnewprice (firstnm) drivetype subsidy tax_sim elast, by(car_option wealth_group)

gen co2_ann=prob_alt*9230/4*et_co2*alt_km/1000
gen tax_total=prob_alt*9230/4*car_tax_alt
gen pv_tax_total=prob_alt*9230/4*pv_tax_alt
gen paid_subsidy=prob_alt*9230/4*subsidy if drivetype==3
gen prob_EV=prob_alt if drivetype==3
gen veh_tax=prob_alt*9230/4*0.04*prhnewprice
replace veh_tax=0 if drivetype==3

gen emb_co2_ann=0
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 + prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_alt*2307.5*alt_km/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gcollapse (sum) co2_ann petrol_tax emb_co2_ann tax_total pv_tax_total paid_subsidy prob_EV veh_tax (firstnm) subsidy tax_sim elast, by(wealth_group)

gen time_fee=3
gen scenario=1
frame put _all, into(optimal`count')

frame change default	
frame drop probs

**Calculate the change in consumer welfare 

frame change default

frame put CS_alt be_alt wealth_group id, into(cs)
frame change cs

gcollapse (firstnm) CS_alt wealth_group (sum) be_alt  , by(id)
gcollapse (sum) CS_alt be_alt , by(wealth_group)

frlink 1:1 wealth_group, frame(cs_init)

frget CS_init be_init, from(cs_init)

gen CS_change=CS_alt-CS_init
gen be_change=be_alt-be_init

replace CS_change=2307.5/5768.5*CS_change

replace be_change=2307.5/5768.5*be_change

frame change optimal`count'
frlink 1:1 wealth_group, frame(cs)
frget CS_change be_change, from(cs)
drop cs 

frame change default
frame drop cs




drop subsidy prob_alt CS_alt be_alt price_alt price_w_alt2 price_w_alt3 price_w_alt4 tax_sim alt_var_cost4 car_tax_alt pv_tax_alt A_K B_K E_K F_K G_K EV_K alt_km elast 


****
*Second case: Elasticity=0, tax_sim=59, Subsidy=5400
****

local count=`count'+1
local y=0
local x=59
local i=5400


**This is the program to calculate the optimal outcomes by quartile: 

gen elast=`y'

gen car_tax_alt=car_tax
replace car_tax_alt=car_tax/0.4 if drivetype==3
replace car_tax_alt=car_tax/0.6 if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax/0.8 if drivetype!=3 & env_cat==2	
	
	gen tax_sim=`x'
	
	merge m:1 tax_sim using "$root/Data/Original/tax_params_ext.dta", 
	drop if _merge==2
	drop _merge
	
	
local r=0.06

gen pv_tax_alt=pv_tax

*for EV
replace pv_tax_alt=car_tax_alt*(1+EV_K)+(car_tax_alt*(1+EV_K))/((1+`r')^(1))+(car_tax_alt*(1+EV_K))/((1+`r')^(2))+(car_tax_alt*(1+EV_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype==3 

*for category A
replace pv_tax_alt=car_tax_alt*(1+A_K)+(car_tax_alt*(1+A_K))/((1+`r')^(1))+(car_tax_alt*(1+A_K))/((1+`r')^(2))+(car_tax_alt*(1+A_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

*for category B
replace pv_tax_alt=car_tax_alt*(1+B_K)+(car_tax_alt*(1+B_K))/((1+`r')^(1))+(car_tax_alt*(1+B_K))/((1+`r')^(2))+(car_tax_alt*(1+B_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*for category E
replace pv_tax_alt=car_tax_alt*(1+E_K)+(car_tax_alt*(1+E_K))/((1+`r')^(1))+(car_tax_alt*(1+E_K))/((1+`r')^(2))+(car_tax_alt*(1+E_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==5
*for category F
replace pv_tax_alt=car_tax_alt*(1+F_K)+(car_tax_alt*(1+F_K))/((1+`r')^(1))+(car_tax_alt*(1+F_K))/((1+`r')^(2))+(car_tax_alt*(1+F_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax_alt*(1+G_K)+(car_tax_alt*(1+G_K))/((1+`r')^(1))+(car_tax_alt*(1+G_K))/((1+`r')^(2))+(car_tax_alt*(1+G_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==7

replace car_tax_alt=car_tax_alt*(1+G_K) if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax_alt*(1+F_K) if drivetype!=3 & env_cat==6
replace car_tax_alt=car_tax_alt*(1+E_K) if drivetype!=3 & env_cat==5
replace car_tax_alt=car_tax_alt*(1+B_K) if drivetype!=3 & env_cat==2
replace car_tax_alt=car_tax_alt*(1+A_K) if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax_alt*(1+EV_K) if drivetype==3 


gen ann_cost=drive_cost_ann+car_tax
gen ann_cost_alt=drive_cost_ann+car_tax_alt


gen alt_km=average_use*(1+elast*((ann_cost_alt-ann_cost)/ann_cost))

replace alt_km=average_use if missing(alt_km)


gen alt_cost_drive=(alt_km/100)*drive_cost
gen alt_drive_cost=alt_cost_drive+alt_cost_drive*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost4=(alt_drive_cost+pv_tax_alt)/1000

replace alt_var_cost4=0 if missing(alt_var_cost4)


drop alt_drive_cost alt_cost_drive ann_cost ann_cost_alt

	
		gen subsidy=`i'
	
	gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-subsidy if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

drop price_alt_tot

sort id car_option

*Add the mata program here: 

mata {
X=. ;
st_view(X,.,"price_alt height weight alt_var_cost4 typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4") ;

var_cost=.
st_view(var_cost,.,"alt_var_cost4")

l_fit=exp(X*beta) ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
fit=mean(l_fit')
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

/*
mata drop prob_all be var_cost 
*/
prob_all=. ;
be=. ;
var_cost=. ; 

st_addvar("double", "prob_alt")
st_store(.,"prob_alt",prob')

st_addvar("double", "be_alt")
st_store(.,"be_alt",be_m')


CS=mean((((-1):/b_price):*log(bigtot))')

/*
mata drop bigtot
*/
bigtot=. ; 

st_addvar("double", "CS_alt")
st_store(.,"CS_alt",CS')

/*
mata drop prob CS be_m
*/
prob=.;
CS=.;
be_m=.;

}


frame put prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt drivetype subsidy prhnewprice tax_sim elast car_option wealth_group, into(probs)

frame change probs

gcollapse (mean) prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt prhnewprice (firstnm) drivetype subsidy tax_sim elast, by(car_option wealth_group)

gen co2_ann=prob_alt*9230/4*et_co2*alt_km/1000
gen tax_total=prob_alt*9230/4*car_tax_alt
gen pv_tax_total=prob_alt*9230/4*pv_tax_alt
gen paid_subsidy=prob_alt*9230/4*subsidy if drivetype==3
gen prob_EV=prob_alt if drivetype==3
gen veh_tax=prob_alt*9230/4*0.04*prhnewprice
replace veh_tax=0 if drivetype==3

gen emb_co2_ann=0
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 + prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_alt*2307.5*alt_km/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gcollapse (sum) co2_ann petrol_tax emb_co2_ann tax_total pv_tax_total paid_subsidy prob_EV veh_tax (firstnm) subsidy tax_sim elast, by(wealth_group)

gen time_fee=3
gen scenario=2
frame put _all, into(optimal`count')

frame change default	
frame drop probs

**Calculate the change in consumer welfare 

frame change default

frame put CS_alt be_alt wealth_group id, into(cs)
frame change cs

gcollapse (firstnm) CS_alt wealth_group (sum) be_alt  , by(id)
gcollapse (sum) CS_alt be_alt , by(wealth_group)

frlink 1:1 wealth_group, frame(cs_init)

frget CS_init be_init, from(cs_init)

gen CS_change=CS_alt-CS_init
gen be_change=be_alt-be_init

replace CS_change=2307.5/5768.5*CS_change

replace be_change=2307.5/5768.5*be_change

frame change optimal`count'
frlink 1:1 wealth_group, frame(cs)
frget CS_change be_change, from(cs)
drop cs 

frame change default
frame drop cs


drop subsidy prob_alt CS_alt be_alt price_alt price_w_alt2 price_w_alt3 price_w_alt4 tax_sim alt_var_cost4 car_tax_alt pv_tax_alt A_K B_K E_K F_K G_K EV_K alt_km elast 


****
*Third case: Elasticity=-0.3, feebate=47, Subsidy=6800
****

local count=`count'+1
local y=-0.3
local x=47
local i=6800


**This is the program to calculate the optimal outcomes by quartile: 

gen elast=`y'

gen car_tax_alt=car_tax
replace car_tax_alt=car_tax/0.4 if drivetype==3
replace car_tax_alt=car_tax/0.6 if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax/0.8 if drivetype!=3 & env_cat==2	
	
	gen tax_sim=`x'
	
	merge m:1 tax_sim using "$root/Data/Original/tax_params_ext.dta", 
	drop if _merge==2
	drop _merge
	
	
local r=0.06

gen pv_tax_alt=pv_tax

*for EV
replace pv_tax_alt=car_tax_alt*(1+EV_K)+(car_tax_alt*(1+EV_K))/((1+`r')^(1))+(car_tax_alt*(1+EV_K))/((1+`r')^(2))+(car_tax_alt*(1+EV_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14))  if drivetype==3 

*for category A
replace pv_tax_alt=car_tax_alt*(1+A_K)+(car_tax_alt*(1+A_K))/((1+`r')^(1))+(car_tax_alt*(1+A_K))/((1+`r')^(2))+(car_tax_alt*(1+A_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

*for category B
replace pv_tax_alt=car_tax_alt*(1+B_K)+(car_tax_alt*(1+B_K))/((1+`r')^(1))+(car_tax_alt*(1+B_K))/((1+`r')^(2))+(car_tax_alt*(1+B_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*for category E
replace pv_tax_alt=car_tax_alt*(1+E_K)+(car_tax_alt*(1+E_K))/((1+`r')^(1))+(car_tax_alt*(1+E_K))/((1+`r')^(2))+(car_tax_alt*(1+E_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==5
*for category F
replace pv_tax_alt=car_tax_alt*(1+F_K)+(car_tax_alt*(1+F_K))/((1+`r')^(1))+(car_tax_alt*(1+F_K))/((1+`r')^(2))+(car_tax_alt*(1+F_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax_alt*(1+G_K)+(car_tax_alt*(1+G_K))/((1+`r')^(1))+(car_tax_alt*(1+G_K))/((1+`r')^(2))+(car_tax_alt*(1+G_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==7

replace car_tax_alt=car_tax_alt*(1+G_K) if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax_alt*(1+F_K) if drivetype!=3 & env_cat==6
replace car_tax_alt=car_tax_alt*(1+E_K) if drivetype!=3 & env_cat==5
replace car_tax_alt=car_tax_alt*(1+B_K) if drivetype!=3 & env_cat==2
replace car_tax_alt=car_tax_alt*(1+A_K) if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax_alt*(1+EV_K) if drivetype==3 


gen ann_cost=drive_cost_ann+car_tax
gen ann_cost_alt=drive_cost_ann+car_tax_alt


gen alt_km=average_use*(1+elast*((ann_cost_alt-ann_cost)/ann_cost))

replace alt_km=average_use if missing(alt_km)


gen alt_cost_drive=(alt_km/100)*drive_cost
gen alt_drive_cost=alt_cost_drive+alt_cost_drive*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost4=(alt_drive_cost+pv_tax_alt)/1000

replace alt_var_cost4=0 if missing(alt_var_cost4)


drop alt_drive_cost alt_cost_drive ann_cost ann_cost_alt

	
		gen subsidy=`i'
	
	gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-subsidy if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

drop price_alt_tot

sort id car_option

*Add the mata program here: 

mata {
X=. ;
st_view(X,.,"price_alt height weight alt_var_cost4 typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4") ;

var_cost=.
st_view(var_cost,.,"alt_var_cost4")

l_fit=exp(X*beta) ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
fit=mean(l_fit')
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

/*
mata drop prob_all be var_cost 
*/
prob_all=. ;
be=. ;
var_cost=. ; 

st_addvar("double", "prob_alt")
st_store(.,"prob_alt",prob')

st_addvar("double", "be_alt")
st_store(.,"be_alt",be_m')


CS=mean((((-1):/b_price):*log(bigtot))')

/*
mata drop bigtot
*/
bigtot=. ; 

st_addvar("double", "CS_alt")
st_store(.,"CS_alt",CS')

/*
mata drop prob CS be_m
*/
prob=.;
CS=.;
be_m=.;

}

frame put prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt drivetype subsidy prhnewprice tax_sim elast car_option wealth_group, into(probs)

frame change probs

gcollapse (mean) prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt prhnewprice (firstnm) drivetype subsidy tax_sim elast, by(car_option wealth_group)

gen co2_ann=prob_alt*9230/4*et_co2*alt_km/1000
gen tax_total=prob_alt*9230/4*car_tax_alt
gen pv_tax_total=prob_alt*9230/4*pv_tax_alt
gen paid_subsidy=prob_alt*9230/4*subsidy if drivetype==3
gen prob_EV=prob_alt if drivetype==3
gen veh_tax=prob_alt*9230/4*0.04*prhnewprice
replace veh_tax=0 if drivetype==3

gen emb_co2_ann=0
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 + prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_alt*2307.5*alt_km/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gcollapse (sum) co2_ann petrol_tax emb_co2_ann tax_total pv_tax_total paid_subsidy prob_EV veh_tax (firstnm) subsidy tax_sim elast, by(wealth_group)

gen time_fee=3
gen scenario=3
frame put _all, into(optimal`count')

frame change default	
frame drop probs

**Calculate the change in consumer welfare 

frame change default

frame put CS_alt be_alt wealth_group id, into(cs)
frame change cs

gcollapse (firstnm) CS_alt wealth_group (sum) be_alt  , by(id)
gcollapse (sum) CS_alt be_alt , by(wealth_group)

frlink 1:1 wealth_group, frame(cs_init)

frget CS_init be_init, from(cs_init)

gen CS_change=CS_alt-CS_init
gen be_change=be_alt-be_init

replace CS_change=2307.5/5768.5*CS_change

replace be_change=2307.5/5768.5*be_change

frame change optimal`count'
frlink 1:1 wealth_group, frame(cs)
frget CS_change be_change, from(cs)
drop cs 

frame change default
frame drop cs




drop subsidy prob_alt CS_alt be_alt price_alt price_w_alt2 price_w_alt3 price_w_alt4 tax_sim alt_var_cost4 car_tax_alt pv_tax_alt A_K B_K E_K F_K G_K EV_K alt_km elast 


****
*Fourth case: Elasticity=-0.3, feebate=30, Subsidy=5600
****

local count=`count'+1
local y=-0.3
local x=30
local i=5600


**This is the program to calculate the optimal outcomes by quartile: 

gen elast=`y'

gen car_tax_alt=car_tax
replace car_tax_alt=car_tax/0.4 if drivetype==3
replace car_tax_alt=car_tax/0.6 if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax/0.8 if drivetype!=3 & env_cat==2	
	
	gen tax_sim=`x'
	
	merge m:1 tax_sim using "$root/Data/Original/tax_params_ext.dta", 
	drop if _merge==2
	drop _merge
	
	
local r=0.06

gen pv_tax_alt=pv_tax

*for EV
replace pv_tax_alt=car_tax_alt*(1+EV_K)+(car_tax_alt*(1+EV_K))/((1+`r')^(1))+(car_tax_alt*(1+EV_K))/((1+`r')^(2))+(car_tax_alt*(1+EV_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype==3 

*for category A
replace pv_tax_alt=car_tax_alt*(1+A_K)+(car_tax_alt*(1+A_K))/((1+`r')^(1))+(car_tax_alt*(1+A_K))/((1+`r')^(2))+(car_tax_alt*(1+A_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

*for category B
replace pv_tax_alt=car_tax_alt*(1+B_K)+(car_tax_alt*(1+B_K))/((1+`r')^(1))+(car_tax_alt*(1+B_K))/((1+`r')^(2))+(car_tax_alt*(1+B_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*for category E
replace pv_tax_alt=car_tax_alt*(1+E_K)+(car_tax_alt*(1+E_K))/((1+`r')^(1))+(car_tax_alt*(1+E_K))/((1+`r')^(2))+(car_tax_alt*(1+E_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==5
*for category F
replace pv_tax_alt=car_tax_alt*(1+F_K)+(car_tax_alt*(1+F_K))/((1+`r')^(1))+(car_tax_alt*(1+F_K))/((1+`r')^(2))+(car_tax_alt*(1+F_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax_alt*(1+G_K)+(car_tax_alt*(1+G_K))/((1+`r')^(1))+(car_tax_alt*(1+G_K))/((1+`r')^(2))+(car_tax_alt*(1+G_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9)) +(car_tax_alt)/((1+`r')^(10)) +(car_tax_alt)/((1+`r')^(11)) +(car_tax_alt)/((1+`r')^(12)) +(car_tax_alt)/((1+`r')^(13)) +(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==7

replace car_tax_alt=car_tax_alt*(1+G_K) if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax_alt*(1+F_K) if drivetype!=3 & env_cat==6
replace car_tax_alt=car_tax_alt*(1+E_K) if drivetype!=3 & env_cat==5
replace car_tax_alt=car_tax_alt*(1+B_K) if drivetype!=3 & env_cat==2
replace car_tax_alt=car_tax_alt*(1+A_K) if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax_alt*(1+EV_K) if drivetype==3 


gen ann_cost=drive_cost_ann+car_tax
gen ann_cost_alt=drive_cost_ann+car_tax_alt


gen alt_km=average_use*(1+elast*((ann_cost_alt-ann_cost)/ann_cost))

replace alt_km=average_use if missing(alt_km)


gen alt_cost_drive=(alt_km/100)*drive_cost
gen alt_drive_cost=alt_cost_drive+alt_cost_drive*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost4=(alt_drive_cost+pv_tax_alt)/1000

replace alt_var_cost4=0 if missing(alt_var_cost4)


drop alt_drive_cost alt_cost_drive ann_cost ann_cost_alt

	
		gen subsidy=`i'
	
	gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-subsidy if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

drop price_alt_tot

sort id car_option

*Add the mata program here: 

mata {
X=. ;
st_view(X,.,"price_alt height weight alt_var_cost4 typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4") ;

var_cost=.
st_view(var_cost,.,"alt_var_cost4")

l_fit=exp(X*beta) ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
fit=mean(l_fit')
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

/*
mata drop prob_all be var_cost 
*/
prob_all=. ;
be=. ;
var_cost=. ; 

st_addvar("double", "prob_alt")
st_store(.,"prob_alt",prob')

st_addvar("double", "be_alt")
st_store(.,"be_alt",be_m')


CS=mean((((-1):/b_price):*log(bigtot))')

/*
mata drop bigtot
*/
bigtot=. ; 

st_addvar("double", "CS_alt")
st_store(.,"CS_alt",CS')

/*
mata drop prob CS be_m
*/
prob=.;
CS=.;
be_m=.;

}


frame put prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt drivetype subsidy prhnewprice tax_sim elast car_option wealth_group, into(probs)

frame change probs

gcollapse (mean) prob_alt et_co2 alt_km et_verbrauch el_verbrauch pv_tax_alt car_tax_alt prhnewprice (firstnm) drivetype subsidy tax_sim elast, by(car_option wealth_group)

gen co2_ann=prob_alt*9230/4*et_co2*alt_km/1000
gen tax_total=prob_alt*9230/4*car_tax_alt
gen pv_tax_total=prob_alt*9230/4*pv_tax_alt
gen paid_subsidy=prob_alt*9230/4*subsidy if drivetype==3
gen prob_EV=prob_alt if drivetype==3
gen veh_tax=prob_alt*9230/4*0.04*prhnewprice
replace veh_tax=0 if drivetype==3

gen emb_co2_ann=0
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_alt*2307.5*alt_km/100*el_verbrauch*0.139 + prob_alt*2307.5*alt_km/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_alt*2307.5*alt_km/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gcollapse (sum) co2_ann petrol_tax emb_co2_ann tax_total pv_tax_total paid_subsidy prob_EV veh_tax (firstnm) subsidy tax_sim elast, by(wealth_group)

gen time_fee=3
gen scenario=4
frame put _all, into(optimal`count')

frame change default	
frame drop probs

**Calculate the change in consumer welfare 

frame change default

frame put CS_alt be_alt wealth_group id, into(cs)
frame change cs

gcollapse (firstnm) CS_alt wealth_group (sum) be_alt  , by(id)
gcollapse (sum) CS_alt be_alt , by(wealth_group)

frlink 1:1 wealth_group, frame(cs_init)

frget CS_init be_init, from(cs_init)

gen CS_change=CS_alt-CS_init
gen be_change=be_alt-be_init

replace CS_change=2307.5/5768.5*CS_change

replace be_change=2307.5/5768.5*be_change

frame change optimal`count'
frlink 1:1 wealth_group, frame(cs)
frget CS_change be_change, from(cs)
drop cs 

frame change default
frame drop cs




drop subsidy prob_alt CS_alt be_alt price_alt price_w_alt2 price_w_alt3 price_w_alt4 tax_sim alt_var_cost4 car_tax_alt pv_tax_alt A_K B_K E_K F_K G_K EV_K alt_km elast 



forvalues i=1/4 {
	frame change optimal`i'
	expand 2 if _n==1, gen(newvar)
	replace wealth_group=0 if newvar==1
	
foreach var of varlist co2_ann emb_co2_ann petrol_tax veh_tax tax_total pv_tax_total paid_subsidy CS_change be_change {
	
	replace `var'=. if wealth_group==0
	egen tot_`var'=sum(`var') 
	replace `var'=tot_`var' if wealth_group==0
	drop tot_`var'	
	
}	
	replace prob_EV=. if wealth_group==0
	egen m_prob=mean(prob_EV) 
	replace prob_EV=m_prob if wealth_group==0
	drop m_prob	

}

frame change optimal1


forvalues ii=2/4 {
	
	frameappend optimal`ii', drop
}

frame rename optimal1 optimal_byquart

**Create the total of subsidies: 
bys subsidy elast tax_sim time_fee: egen tot_subs=total(paid_subsidy)

*Calculate the share of subsidy payments: 
gen subs_share=100*paid_subsidy/tot_subs

*Calculate the share of levy payments:
bys subsidy elast tax_sim time_fee: egen tot_tax=total(pv_tax_total)

gen tax_share=100*pv_tax_total/tot_tax

sort elast subsidy tax_sim time_fee wealth_group
save "$root/Results/estimation_results/opt_pol_byquart.dta", replace

